#!/usr/bin/env python
# -*- coding: UTF-8 -*-
# Edit By WangWenHao, CMIOT, China Mobile

import sys
import time
import pandas as pd
import numpy as np
from datetime import datetime
import matplotlib.pylab as plt

#时间格式转为分钟
def time2min(time_str):
    h,m = time_str.strip().split(":")
    return int(h) * 60 + int(m)

#分钟转为时间格式
def min2time(min_cnt):
	h, m = divmod(min_cnt, 60)
	str = "%d:%d" % (h,m)
	return str

#自定义的画图
def myplot(data, title, xlab, ylab):
	plt.rcParams['font.sans-serif']=['SimHei']
	plt.rcParams['axes.unicode_minus'] = False
	plt.plot(data)
	plt.title(title)
	plt.xlabel(xlab)
	plt.ylabel(ylab)
	plt.show()

############################################################################################
########       quote: https://www.freesion.com/article/326951832/                  #########
#平稳性检验
def test_stationarity(timeseries):
    # Determing rolling statistics
    from statsmodels.tsa.stattools import adfuller
    rolmean = timeseries.rolling(window=12, center=False).mean()
    rolstd =  timeseries.rolling(window=12, center=False).std()
    
    # Plot rolling statistics:
    orig = plt.plot(timeseries, color='blue', label='Original')
    mean = plt.plot(rolmean, color='red', label='Rolling Mean')
    std = plt.plot(rolstd, color='black', label = 'Rolling Std')
    plt.legend(loc='best')
    plt.title('Rolling Mean & Standard Deviation')
    plt.show(block=False)
    
    # Perform Dickey-Fuller test:
    print('Results of Dickey-Fuller Test:')
    dftest = adfuller(timeseries, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
    for key,value in dftest[4].items():
        dfoutput['Critical Value ({})'.format(key)] = value
    print(dfoutput)
############################################################################################

# 画acf、pacf图
def acf_draw(ts):
	from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
	plot_acf(ts, lags=30)    #man lags 30
	plt.show()
	plot_pacf(ts, lags=30)   #max lags 30
	plt.show()

#adf单位根检验
def adf_test(ts):
	from statsmodels.tsa.stattools import adfuller
	print(adfuller(ts))

#kpss单位根检验
def kpss_test(ts):
	from statsmodels.tsa.stattools import kpss
	print(kpss(ts))

#Arima模型拟合
def arima_fit(ts, p, d, q, show):
	from statsmodels.tsa.arima_model import ARIMA
	arima_model = ARIMA(ts, (p,d,q)).fit()
	print("arima_model(%d,%d,%d) AIC = %d" %(p ,d, q, arima_model.aic))
	if show == True:
		plt.plot(ts)
		plt.plot(arima_model.fittedvalues, color='red')
		plt.title('RSS: %.4f'% sum((arima_model.fittedvalues-ts)**2))
		plt.show()
	else:
		pass
	return arima_model

#源数据文件读取
def prepare():
	global df
	df = pd.read_csv('../newdata_1.csv', sep = ";")
	df.set_index('newindex', inplace = True)
	#print(df)

#生成睡眠时长时间序列
def MakeTs_bedtime():
	#global ts_bedtime
	df.index = pd.to_datetime(df.index)  # turn index to ts
	ts_bedtime = df['Time in bed'].astype('float64') 
	#print("Head of ts_bedtime:") 
	#print(ts_bedtime.head())
	return ts_bedtime

#生成睡眠质量时间序列
def MakeTs_quality():
	#global ts_quality
	df.index = pd.to_datetime(df.index)  # turn index to ts
	ts_quality = df['Sleep quality'].astype('float64')
	#print("Head of ts_quality:")
	#print(ts_quality.head())
	return ts_quality

#1阶差分
def diff(ts):
	return ts.diff(1)

#指数平滑
def log_trans(ts):
	return np.log(ts)

#移动平滑
def wei_rol_mean(ts, size):
    # 对size个数据进行加权移动平均
    rol_weighted_mean = pd.ewma(ts, span=size)
    return rol_weighted_mean	

#相关性分析
def corr_analyze():
	df = pd.read_csv('../newdata_1.csv', sep = ";")
	df.drop(['newindex'],axis=1,inplace=True)
	df.drop(['End'],axis=1,inplace=True)
	df.drop(['Heart rate'],axis=1,inplace=True)
	mins = []
	for i in range(0, len(df)):
		string = df['Start'][i]
		hour = int(string[11])*10 + int(string[12])
		miniute = int(string[14])*10 + int(string[15])
		mins.append(60*hour+ miniute)
	df.drop(['Start'],axis=1,inplace=True)
	df['Start_mins'] = (mins)	
	print(df.corr())




if __name__ == '__main__':
	# read the newdata_0 as df

	#corr_analyze()

	#sys.exit(0) #STOP HERE

	print("Forecast prepare...")
	prepare()
	ts_bedtime = MakeTs_bedtime()
	ts_quality = MakeTs_quality()
	diff_ts_qua = diff(ts_quality)
	log_ts_qua = log_trans(ts_quality)
	mean_ts_qua = wei_rol_mean(ts_quality, 7)
	ts_train = wei_rol_mean(ts_quality, 7)
	#ts_for_train = mean_ts_qua[0:830]
	#myplot(ts_for_train, '睡眠质量', '日期', '百分比')
	#print(mean_ts_qua)
	#diff_ts_qua[0] = 0.0
	#print(diff_ts_qua)
	#myplot(ts_bedtime, '睡眠时间', '日期', '分钟')
	#myplot(ts_quality, '睡眠指令', '日期', '百分比')
	#test_stationarity(ts_bedtime)
	#test_stationarity(ts_quality)
	#acf_draw(ts_quality)
	#myplot(mean_ts_qua, '睡眠质量', '日期', '百分比')
	#adf_test(mean_ts_qua)
	#kpss_test(mean_ts_qua)
	#acf_draw(mean_ts_qua)
	#adf_test(ts_bedtime)
	#kpss_test(ts_bedtime)
	#acf_draw(ts_bedtime)
	#sys.exit(0)   #STOP HERE
	##########  对睡眠质量序列进行线性拟合 并评估精度   #########
	if 0:
		arima_quality = arima_fit(ts_train, 1, 0, 0, False)
		arima_quality = arima_fit(ts_train, 0, 0, 1, False)
		arima_quality = arima_fit(ts_train, 1, 0, 1, False)
		arima_quality = arima_fit(ts_train, 2, 0, 1, False)
		arima_quality = arima_fit(ts_train, 2, 0, 0, False)
		arima_quality = arima_fit(ts_train, 3, 0, 0, False)
		arima_quality = arima_fit(ts_train, 4, 0, 0, False)


	#sys.exit(0)   #STOP HERE
	fore_ts = []
	for i in range(0, 30):
		arima_quality = arima_fit(ts_train[0:830+i], 1, 0, 0, False)
		fore_data = arima_quality.forecast(steps = 30)[0]
		print(ts_train[831])
		print(fore_data[0])
		ts_train[831+i] = fore_data[0]

	plt.figure(facecolor='white')
	ts_train[830:860].plot(color='red', label='Predict')
	mean_ts_qua[800:860].plot(color='blue', label='Original')
	plt.legend(loc='best')
	#计算 RMSE 评估拟合准确度
	plt.title('RMSE: %.4f'% np.sqrt(sum((ts_train[830:860]-mean_ts_qua[830:860])**2)/mean_ts_qua[830:860].size))
	plt.show()
	#sys.exit(0)   #STOP HERE

	rol_sum = ts_train.rolling(window=6).sum()
	rol_recover = ts_train*7- rol_sum.shift(1)

	plt.figure(facecolor='white')
	rol_recover[830:860].plot(color='red', label='Predict')
	ts_quality[800:860].plot(color='blue', label='Original')
	plt.legend(loc='best')
	#计算 RMSE 评估拟合准确度
	plt.title('RMSE: %.4f'% np.sqrt(sum((rol_recover[830:860]-ts_quality[830:860])**2)/ts_quality[830:860].size))
	plt.show()


##########